home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / rk4.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  960 b   |  43 lines
  1. /*
  2. ### 4th orderRunge-Kutta integrator one-stepper ###
  3.  
  4. Note:     slightly modified version of 4-th order Runge-Kutta integrator "rk4"
  5.     in "Numerical Recipes"
  6.     c convention of [0,n-1] is used
  7. */
  8.  
  9. void
  10. rk4(y,dydx,n,x,h,yout)
  11. double y[],dydx[],x,h,yout[];
  12. int n;
  13. {
  14.     int i;
  15.     double xh,hh,h6,*dym,*dyt,*yt,*dvector();
  16.     void free_dvector();
  17.     extern int var_dim,model;
  18.     extern double *param;
  19.     extern int (*f_p)();
  20.  
  21.     dym = dvector(0,n-1);
  22.     dyt = dvector(0,n-1);
  23.     yt = dvector(0,n-1);
  24.     hh = h * 0.5;
  25.     h6 = h / 6.0;
  26.     xh = x + hh;
  27.     for(i=0;i<n;i++) yt[i] = y[i] + hh * dydx[i];
  28.     (int) f_p(dyt,0,yt,param,xh,var_dim);
  29.     for(i=0;i<n;i++) yt[i] = y[i] + hh * dyt[i];
  30.     (int) f_p(dym,0,yt,param,xh,var_dim);
  31.     for(i=0;i<n;i++){
  32.         yt[i] = y[i] + h * dym[i];
  33.         dym[i] += dyt[i];
  34.     }
  35.     (int) f_p(dyt,0,yt,param,x+h,var_dim);
  36.     for(i=0;i<n;i++){
  37.         yout[i] = y[i] + h6 * (dydx[i] + dyt[i] + 2.0 * dym[i]);
  38.     }
  39.     free_dvector(yt,0,n-1);
  40.     free_dvector(dyt,0,n-1);
  41.     free_dvector(dym,0,n-1);
  42. }
  43.